In [1]:
%matplotlib inline
import numpy as np
import scipy.stats as stat
from matplotlib import pyplot as plt
Generating coupon probabilities can be accomplished in much the same manner as with the normal probabilities. It is actually a little easier because the exponential distribution is defined by only one parameter ($\lambda$, the rate parameter). The mean of the exponential distribution is $\frac{1}{\lambda}$, and we know that the mean coupon probability must be $\frac{1}{n}$, therefore we must set $\lambda=n$. We can generate the probabilities in much the same way as before:
In [2]:
n = 20 #number of coupons
scale = 1/n #scipy uses the scale parameter instead of lambda. Scale and lambda are reciprocals of each other.
x = np.arange(n)+0.5 #arange goes from 0 to n-1, and I want it to go from 1 to n
p_x = stat.expon.ppf(x/n, loc=0, scale=scale)
print('unfilled probability: ', 1-np.sum(p_x))
p_x = p_x/np.sum(p_x)
print('mean prob = %1.16f' % np.mean(p_x), ', sum of all probs = %1.16f' % np.sum(p_x) )
plt.bar(x, p_x, align='center')
plt.axis([0,n,0,4/n])
Out[2]:
In generating the probabilities for the exponential distribution we had to re-scale the probabilities so they would sum to 1. This is because of the "discretization" of the distribution leaves a little bit of the probabilities unfilled. It's not a lot (1.7% missing when n=20), but it needs to be accounted for so we'll re-scale the distribution. I'm not sure of another way to generate exponential probabilities without having to re-scale afterwards.
As the population size increases, we do not approach a uniform distribution, unlike the normal distribution. The shape remains "exponential," low numbered coupons remain rare compared to high numbered coupons.
In [3]:
n = 2000 #number of coupons
scale = 1/n #scipy uses the scale parameter instead of lambda. Scale and lambda are reciprocals of each other.
x = np.arange(n)+0.5 #arange goes from 0 to n-1, and I want it to go from 1 to n
p_x = stat.expon.ppf(x/n, loc=0, scale=scale)
print('unfilled probability: ', 1-np.sum(p_x))
p_x = p_x/np.sum(p_x)
print('mean prob = %1.16f' % np.mean(p_x), ', sum of all probs = %1.16f' % np.sum(p_x) )
plt.bar(x, p_x, align='center')
plt.axis([0,n,0,4/n])
Out[3]: